home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_4 / a4_3.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  3.5 KB  |  146 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 4.3 (Lagrange Approximation).
  9. % Section    4.3,    Lagrange Approximation, Page 224
  10. echo on; clc; hold off; clear
  11. % Investigation of Lagrange polynomial approximations.
  12.  
  13. % n+1  are points needed to construct  Pn(x).
  14.  
  15. % The abscissas are stored in  X.
  16.  
  17. % The ordinates are stored in  Y.
  18.  
  19. % The points are counted  k=1,2,...,n+1.
  20.  
  21. % The Lagrange coefficient polynomial L   (x) 
  22. %                                      n,k   
  23.  
  24. % is stored in row  k  of the matrix  L(k,j)
  25.  
  26. % Remark. lagran.m is used for Algorithm 4.3
  27.  
  28. pause % Press any key to continue.
  29.  
  30. clc;
  31. % Lagrange polynomial approximation Pn(x) of f(x) over [a,b].
  32.  
  33. % Place the degree of approximation in  n.
  34.  
  35. % Place the left endpoint in  a.
  36.  
  37. % Place the right endpoint in b.
  38.  
  39. % Place the 'string expression' for f(x) in  fun.
  40.  
  41. n = 3;
  42. a = 0;
  43. b = 1.2;
  44. fun = 'cos(x)';
  45.  
  46. pause % Press any key to form the coefficient polynomials.
  47.  
  48. clc;
  49. % The Lagrange interpolating polynomial is being constructed.
  50.  
  51. h = (b-a)/n;
  52.  
  53. X = a:h:b;
  54.  
  55. x = X;
  56. Y = eval(fun);
  57.  
  58. C = lagran(X,Y);
  59.  
  60. pause % Press any key to continue.
  61.  
  62. clc;
  63. format short;
  64. Mx1 = 'Construction of a Lagrange approximation polynomial.';
  65. Mx2 = 'The abscissas are:';
  66. Mx3 = 'The ordinates are:';
  67. clc,echo off,diary output,...
  68. disp(''),disp(Mx1),disp(''),disp(Mx2),disp(X),disp(Mx3),disp(Y),...
  69. diary off, echo on
  70.  
  71. pause  % Press any key to continue.
  72.  
  73. clc;
  74. format short;
  75. Mx1 = 'Construction of a Lagrange approximation polynomial.';
  76. Mx2 = 'The abscissas are:';
  77. Mx3 = 'The ordinates are:';
  78. clc,echo off,diary output,...
  79. disp(''),disp(Mx1),disp(''),disp(Mx2),disp(X),disp(Mx3),disp(Y),...
  80. diary off, echo on
  81.  
  82. pause  % Press any key graph  f(x)  and  Pn(x).
  83.  
  84. clc; clg;
  85. a = min(X);
  86. b = max(X);
  87. h = (b-a)/150;
  88. X1 = a:h:b;
  89. Y1 = polyval(C,X1);
  90. x = X1;
  91. Y2 = eval(fun);
  92. plot(X,Y,'or',X1,Y1,'-r',X1,Y2,'-g');...
  93. hold on;...
  94. plot([a b],[0 0],'b',[0 0],[-10000 10000],'b');...
  95. xlabel('x');...
  96. ylabel('y');...
  97. Mx1 = ['Comparison of ',fun,' and P'];...
  98. Mx2 = [Mx1,num2str(n),'(x).'];...
  99. title(Mx2);...
  100. grid;...
  101. axis;...
  102. hold off;...
  103. shg; pause    % Press any key to continue.
  104.  
  105. Mx1='The Lagrange polynomial has been rearranged in ordinary polynomial form.';
  106. Mx2='This ordinary polynomial looks like:';
  107. Mx3='Pn(x) = c(1)x^n + c(2)x^(n-1) + ... + c(n)x + c(n+1)';
  108. Mx4 = 'The degree is  n = ';
  109. Mx5 = ', and the coefficients list  C  is:';
  110. clc,echo off, diary output,...
  111. disp(''),disp(Mx1),disp(''),disp(Mx2),disp(''),...
  112. disp(Mx3),disp(''),disp([Mx4,num2str(n),Mx5]),disp(''),...
  113. for i=1:5:n+1, disp(C([i:min(i+4,n+1)])); end,...
  114. diary off, echo on
  115.  
  116. pause  % Press any key to view f(x) - Pn(x).
  117.  
  118. clc; clg;
  119. n1 = length(X);
  120. Z = zeros(1,n1);
  121. plot(X,Z,'or',X1,Y2-Y1,'-r');...
  122. hold on;...
  123. plot([a b],[0 0],'b',[0 0],[-10000 10000],'b');...
  124. xlabel('x');...
  125. ylabel('y');...
  126. Mx1 = ['The error: ',fun,' - P'];...
  127. Mx2 = [Mx1,num2str(n),'(x).'];...
  128. title(Mx2);...
  129. grid;...
  130. hold off;...
  131. shg; pause    % Press any key to continue.
  132.  
  133. pause    % Press any key for a list of numerical computations.
  134.  
  135. clc;
  136. X = a:0.1:b;
  137. x = X;
  138. Y = eval(fun);
  139. P = polyval(C,X);, format long;
  140. points = [X;Y;P;Y-P];
  141. Mx1=['Lagrange polynomial approximation of f(x) = ',fun];
  142. Mx2='     x(k)               f(x(k))            Pn(x(k))           error';
  143. clc,echo off,diary output,...
  144. disp(''),disp(Mx1),disp(''),disp(Mx2),disp(points'),...
  145. diary off,echo on
  146.